Semiclassical description of spin ladders 

D. Senechal 

Centre de Recherche en Physique du Soltde et Departement de Physique, 
Universite de Sherbrooke, Sherbrooke, Quebec, Canada J1K 2R1. 
(June 1995) 

The Heisenberg spin ladder is studied in the semiclassical limit, via a mapping to the nonlinear a model. Different 
treatments are needed if the inter-chain coupling K is small, intermediate or large. For intermediate coupling a 
single nonlinear a model is used for the ladder. Its predicts a spin gap for all nonzero values of K if the sum 
s + s of the spins of the two chains is an integer, and no gap otherwise. For small K, a better treatment proceeds 
by coupling two nonlinear sigma models, one for each chain. For integer s — s, the saddle-point approximation 
predicts a sharp drop in the gap as K increases from zero. A Monte-Carlo simulation of a spin 1 ladder is 
presented which supports the analytical results. 



I. INTRODUCTION 

Interest in one-dimensional quantum antiferromagnets 
has been great ever since Haldane conjecturedEl that 
integer-spin chains have a gap in their excitation spec- 
trum while half-integer spin chains do not. More recently, 
coupled spin chains have aroused some interest, in partic- 
ular the so-called Heisenberg ladder, whose Hamiltonian 
may be written as follows: 

H = J £ {Si • Sj+i + S s: • S i+1 } + K ]T S, • S, (1) 

i i 

wherein Si and Si are the spin operators at site i for 
chains A and B respectively. The relevant parameters are 
the spins s and s of each chain, and the coupling ratio 
K / J . Since this system is one-dimensional and is char- 
acterized by a continuous order parameter, one does not 
expect any long range order for any value of the couplings 
J and K. The important question is rather how the spin 
excitation gap A develops or varies as the inter-chain cou- 
pling K is turned on. The case of two spin-^ chains is the 
one that has received attention recently, and has been 
studied with various techniques: BosonisationQ, exact 
diagonalizationaa, Quantum Monte-CarloS and density- 
matrix renormalizatioroQ. The prevailing conclusion is 
that a gap A appears for any nonzero inter-chain cou- 
pling K. This gap increases with K in an almost linear 
fashion. 

In this work we will study the Heisenberg ladder with 
semi-classical techniques, i.e., by using an approximate 
mapping from the Hamiltonian to the nonlinear a model. 
As long as K and J are not too different, one may assume 
short range antiferromagnetic (AF) order along and ac- 
cross the chains, and work out a direct mapping with the 
nonlinear a model. This will be described in section II. 
However, the validity of such a mapping is questionable 
when K/J is too small (or too large). If K <§C J one 
should rather consider two coupled nonlinear a models, 
one for each chain. For spin-^ chains this analysis is 



complicated by the existence of topological terms, and 
we shall carry it only for integer spin; this is done in sec- 
tion III. In section IV we will compare these semi-classical 
results with the outcome of a quantum Monte-Carlo sim- 
ulation of the spin-1 Heisenberg ladder. 

II. INTERMEDIATE COUPLING ANALYSIS 

A. Derivation of the a model 

In this section we will show how the AF Heisenberg 
ladder defined by the Hamiltonian (Q) may be described, 
in the continuum and semi-classical limit, by the one- 
dimensional nonlinear a model with Lagrangian density 

Ca = Tg{\ 1 [dtm)2 ~ V[dxm)2 ) (2) 

wherein the field m(x,t) is a unit vector. The coupling 
g and the velocity v obtained semiclassically depend on 
the spins s and s, and on the inter-chain and intra-chain 
couplings K and J. If the two chains are alike (s = s), 
these parameters are 

v = 2asJy/l + K/2J g= -^/l + K/2J (3) 

s 

If s 7^ s, the corresponding expressions are more compli- 
cated. 

In addition to the a model Lagrangian (|2|), there also 
appears a topological term 

Q 

C to P = • (9 t m x d x m) (4) 

with = 2tt{s — s). One easily shows that, with suitable 
boundary conditions at infinity, the action obtained by 
integrating £ top is an integer times 0. It follows that 
this topological term has no effect if s — s is an integer. 
On the other hand, if s — s is a half-integer, i.e., if one 
of the chains is made of half-integer spins and the other 
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of integer spins, then the topological term appears as it 
does in the semiclassical analysis of a single half-integer 
spin chain, from which we conjecture, following Haldane, 
that such a system does not have a gap. 

Let us now go through this semiclassical derivation. 
We will proceed as in Ref. |ll|. The idea is to first de- 
scribe the dynamics of a single spin in the Lagrangian 
formalism, with the help of spin coherent states. For de- 
tails_Dn this technique the reader is referred to Fradkin's 
textile or to the review by Manousakisliia. The dynamics 
of a single spin S without interaction is described by a 
unit vector n such that S = sn and the corresponding 
action, named after Wess and Zumino (WZ), is nonlocal: 

Swz = I dt C kin = s dr dtn- (d T n x d t n) 
Jo Jo Jo 

(5) 

Here the time t runs from to some finite period T 
and t is an auxiliary coordinate introduced in order to 
parametrize, along with t, the spherical cap delimited 
by the curve n(t) from t = to t = T. What is actually 
needed is the variation of this action upon a small change 
Sn: 

SS WZ — s / dt Sn ■ (n x d t n) (6) 
Jo 

Each spin on the ladder may be described by a unit 
vector n like above. We expect short range AF order, 
with a 'magnetic cell' containing four spins. If we call r 
the position of a reference spin in each magnetic cell, the 
positions r' of each of the four spins within the cell may 
be written as 

r' = r + ax + /3y (7) 

where x and y are the lattice vectors (respectively along 
and across the ladder) and where a and (3 each run from 
to 1. The short range AF order motivates the intro- 
duction of a unit vector n which varies slowly in space: 

S(r) = sn(r) 

S(r + x) = — sn(r + x) 

S(r + y) = -sn(r + y) 

S(r + x + y) = sn(r + x + y) (8) 

The four spins of the unit cell may be described differ- 
ently, at our convenience, by four different vectors, which 
we choose as follows (in terms of the unit vector defined 
above): 

n(r) = m(r) + a [ l i(r) + l 10 (r) + l u (r)] 
n(r + x) = m(r) + a [ l i(r) - li (r) - l u (r)] 
n(r + y) = m(r) + a [-loi(r) + li (r) - l u (r)] 
n(r + x + y) = m(r) + a [-loi(r) - lio(r) + l u (r)] (9) 



The mean value of n within a magnetic cell is m. In 
addition, we have introduced three deviation field 1 10 , 
loi and In which describe deviations from strict AF or- 
der. We have included a factor of the lattice spacing a 
in front of these deviation fields in order to emphasize 
our assumption that the deviations from the short-range 
AF order are small, and to control this approximation in 
the same way as the continuum limit. The four fields m 
and ly all depend on the cell position r = ix only, since 
they are defined for a cell as a whole, and are entirely 
equivalent to the specification of the four spins S(r') of 
that cell. They are not independent, but must obey the 
constraints imposed by the relation n 2 = 1. Thus, we 
have 4x2 = 8 degrees of freedom per unit cell. 

So far the treatment has been exact. Now we will ex- 
press the kinetic term Skin (the sum of the Wess-Zumino 
actions for each of the spins) and the potential term (the 
Heisenberg Hamiltonian) in terms of these new variables, 
to lowest nontrivial order in a. We will then integrate 
out (in the functional sense) the deviation fields to end 
up with a theory defined only in terms of the staggered 
magnetization m. 

Let us start with the kinetic term. In principle this 
term is nonlocal, as expressed in Eq. (pi). The local AF 
order allows us to write it in local form (to first order in a) 
since the WZ action is odd under the reversal of spin and 
Eq. (||) may be applied between two neighboring sites. If 
we define the difference 

S x n(r') = n(r' + x) - n(r') (10) 

(likewise for S y n) one may write the kinetic term, to first 
order in a, as 

Skin = ^2 [sS x n(r) + sS x n(r + y)] • (m x <9 t m) (11) 

r 

Since 5 x n(r) = — 2a(Zio + iii) and 8 x n(r + y) = — 2a(iio — 
in), this reduces, after replacing the sum by an integral, 
to the following: 

Skin = J dxdt [(s + s)ln + (s - s)l w ] ■ (m X <9 t m) 

(12) 

The Heisenberg Hamiltonian may be separated into 
three parts: H = V x + V x + V y where, up to an additive 
constant, 

r 

v* = I s2 3 E { [^ n ( r + y)] 2 + t^ n ( r + y + *)1 2 } ( 13b ) 

r 

V y = issA^l^ntr)] 2 + [<S„n(r + x)] 2 } (13c) 



Expressing the above in terms of the fields defined in 




we perform a first order Taylor expansion on m, while we 
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essentially neglect the spatial variations of the deviation 
fields. For instance, we have 



m • L 







(20) 



S x ii(r + y + x) = 2ad x m + 2a(l w - In) 



(14) 



Converting the sums of Eq. ( |l3a|) into integrals, we find 
the following Hamiltonian density 



H = 2Jas 2 



2Jas 



7.2 



(d x m) 2 + (li + In) • (1 10 + In + d x m) 
(d x m) 2 + (ho - In) • (lie - In + d x m) 



2Kass(l 2 Q1 +l 2 11 ) 



(15) 



The complete Lagrangian density is of course C = Ckin — 

n. 

The next step is the functional integration of the devia- 
tion fields. Since these fields occur at most in a quadratic 
fashion, they may be integrated simply by solving the 
classical equations of motion and substituting the solu- 
tions back into the Lagrangian density. Taking variations 
of C with respect to lio, loi and In, one finds, after some 
algebra, 



hi = -7—? — — (m x o t m) 



AaJ ss 



lw = -l dxm+imXdtm) ij^V ' 2J s~s 
Substituting into the expression for £, one finds 



K (s + s) 5 



(16a) 
(16b) 

(16c) 



C = -(s-s)m ■ (d t m x 8 x m) + A(d t m) 2 - B(d x m) 2 



where the constants A and B are 

1 \ + {K/2J){s-S) 2 /AsS 



A 



B 



4a J 1 + {K/2J){s 2 + P)/2ss 
l -aJ{s 2 + S 2 ) 



(17) 

(18a) 
(18b) 



In this form the constraints are entirely compatible with 
the equations ([!(]) and we may forget our misgivings. The 
constraint m 2 = 1 is part of the definition of the model 
(|J). The characteristic velocity of Eq. (||) is then v — 
y/B/A and the coupling is g = 1/2 Av. The coefficient of 
the topological term is indeed 9 = 2tt(s—s), as announced 
above. It is easily checked that the correct expressions 
for v and g in the case s — s coincide with Eq. ([}]) . This 
concludes the semi-classical derivation of the nonlinear 
model for the Heisenberg ladder. 



B. Estimating the gap 

The 0(3)/0(2) nonlinear a model - as the model de- 
fined in (Q) is precisely known - lends itself to an estimate 
of the spin excitation gap in the saddle-point approxima- 
tion. This has been described in the literature (cf. the 
reviews in Refs. |l4]and[l^) but will be summarized briefly 
here. 

If it were not for the constraint m 2 = 1, the model 
defined by the Lagrangian (^) would simply describe 
free, massless excitations. A classical analysis of the 
model (incorporating the constraint) reveals a sponta- 
neous breaking of the internal rotation symmetry to a 
ground state in which the field m is uniform. However, 
this ordered state is destroyed by quantum fluctuations 
in dimension 1. This is revealed by the following saddle- 
point analysis: First, one implements the unit-vector con- 
straint by introducing a Lagrange multiplier a(x, t) at ev- 
ery space-time point and by adding the term <r(m 2 — 1) to 
the Lagrangian density. After rescaling the fields and go- 
ing to imaginary time r = —it, the complete Lagrangian 
may be written as 



C E = l -{d m ) 2 + l -a{ m 2 -l/g) 



(21) 



Of course, a simple substitution of Eqs ( |16[ ) into C 
may seem dishonest since, as we indicated earlier, the 
deviations fields are not independent, but constrained by 
the fixed length of each spin. Expressed in terms of m 



Here (9m) 2 stands for (d T m) 2 + (d x m) 2 . We have set the 
velocity v equal to one and will restore it later by dimen- 
sional analysis. The field m is now unconstrained, and its 
functional integration yields the effective potential V(a) 
for the multiplier field a. In the saddle-point approxi- 
mation, the fluctuations of a are neglected (cr is then a 
constant) and the derivative of the effective potential is 



and ljj, these constraints are 



m 2 + a 2 (l 2 1 +l 2 + l 2 1 ) = l 
m • l i + alio • In = 
m • lio + al i ■ In = 
m • In + alio • loi = 



V'(a) = 



(19a) 
(19b) 
(19c) 
(19d) 



1 
1 

2^ 



3 
2 

a 

8tt 



dfc duj 



1 



27T 27T UJ 2 +k 2 + <T 



ln(cr/A 2 



(22) 



Since the classical equations ( |16| ) are valid only at low- 
est order in a, it suffices to apply the above constraint 
equations at that order: 



wherein A is a momentum cutoff, such that Aa ~ 1. 
The saddle-point is thus fixed by the vanishing of the 
above, which occurs at <To = A 2 exp — (An/3g). A nonzero 
value of the saddle-point means that the field a should 
be redefined with respect to the position of the saddle, 
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as a = a + erg, wherein a is a field fluctuating about zero 
and (To is a constant multiplying |m 2 in the quantum 
effective action of the nonlinear er model. Thus, <tq has 
the interpretation of a quantum fluctuation-induced mass 
squared for the field m. At this level of approximation 
we consider no other correction to the effective action: 
in the language of the 1/N expansion, we stay at lowest 
order. Restoring the characterstic velocity, the excitation 
gap may be expressed as 



A = Auexp-(27r/3g) 



(23) 



Note that this expression is 'nonpcrturbative', in the 
sense that it admits no series expansion in powers of g 
about 3 = 0. Since g ~ 1/s is the natural expansion pa- 
rameter of spin-wave theory, it is quite understandable 
that the latter cannot account for the existence of the 
Haldane gap. 



A/A, 



K/J 



FIG. 1. Reduced gap A/Aq as a function of inter-chain 
coupling K/J, as obtained in the saddle-point approximation 
of the non-linear sigma model. The value of Ao comes from 
the saddle-point approximation in the single-chain nonlinear 
sigma model. 

If we now substitute the expression (|3|) for v and g, we 
find 



A/ J = 2(Aa)s y / l + K/2Jexp- 



2-ks 



3^1 + K/2J 



(24) 



This formula indicates a growth of A as a function of K, 
due to the variation of both g and v. Since the model is 
defined on the continuum, an arbitrary lattice spacing a 
has been introduced and the precise value of the gap can- 
not be known. However, it is hoped that its dependence 
on various parameters may be well approximated by this 
method. The /^-dependence of the gap A obtained in 
this fashion is illustrated on Fig. [I]. Up to fairly large 
values of K, the dependence is almost linear. The small 
K behavior is grossly wrong, as will be discussed below. 



C. Extensions 

The calculation of part IIA may also be applied to 
the case of a ferromagnetic inter-chain coupling (K < 
0). We simply have to redefine the spins of the second 
chain by a sign: Si — > — S%. Of course Eq. (@) has to 
be modified (the last two r.h.s. change sign) but Eq. (||) 
stays the same. We find that the Hamiltonian density 
( |l5| ) is unchanged, but the kinetic term is modified to 



Ckin = 2slio • (m x d t m) 



(25) 



(we set s = s for simplicity). Integrating out the devia- 
tion fields, we find in this case that all traces of K disap- 
pear! We again find the nonlinear a model Lagrangian 



but this time with the parameters 



2Jas 



1 

9=- 

s 



(K < 0) 



(26) 



If we compare with the corresponding parameters for a 
single isolated chain (v — 2 J as and g = 2/s), we see 
that the ladder behaves as if it were a single chain of 
effective spin 2s, albeit with a characteristic velocity half 
as large, which implies a gap half as large as that of the 
single chain. In other words, this approach implies that 
the spin-i ladder with ferromagnetic inter-chain coupling 
should behave like a single spin-1 chain, but with half the 
usual Haldane gap: A as 0.205 J, for all values of K. Of 
course, like everything in this section, this result is not to 
be trusted when K <C J ■ However, it has been confirmed, 
by density-matrix renormalization group calculationsE£l 
for reasonable values of K. 

Another possible extension of the above analysis is a 
variant of the Kondo necklace, wherein the coupling is 
taken to be Heisenberg-like. The corresponding Hamil- 
tonian is obtained from (|I|) by deleting the • S-j+i term. 
In this system we also expect short range AF order and 
we may use the same decomposition (JoJ) . Here the Hamil- 
tonian density will be simpler (the second line of Eq. jl5| ) 
will disappear) with the end result v = V JK as and 
g={\/a)y/KfJ. 



III. WEAK COUPLING ANALYSIS 

The results of the previous section are no longer appli- 
cable if K <C J ■ Indeed, some of the heuristic assump- 
tions made in the derivation of the nonlinear a model 
are violated: a small inter-chain coupling implies that 
inter-chain AF correlations are weaker and the deviation 
field aloi is not necessarily small. We may guess the 
rough effect this will have on the gap by the following 
crude argument. Assuming that the nonlinear a model 
is still applicable, the fact that aloi is not small alters the 
constraint on the staggered magnetization. After rescal- 
ing the fields as in Eq. (|l]), this constraint would be 
m 2 = l/g e ff < 1/.9- The 'effective coupling' g e ff is 
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greater than q, and accordingly the gap is larger, as im- 
plied by Eq. (£3|). We thus would expect the actual gap 
to be greater than what is predicted by Eq. (24), as K 
approaches zero. In fact, the gap obtained for a single 
spin chain in the saddle-point approximation of the non- 
linear sigma model is A/ J = 2s(Aa)e~ 7rs / 3 , whereas the 
K -> limit of Eq. (§|) is smaller by a factor e™/ 3 . Of 
course, this holds for integer spin only. For half-integer 
spin chains, the gap should disappear as K — > 0, accord- 
ing to the Haldane conjecture. 

Since the deviation aloi is not necessarily small at weak 
coupling, the constraints ( |l9| ) can no longer be simplified 
and the integration of the deviation fields is no longer 
straightforward. Then we might as well treat the two 
chains separately in the semi-classical approximation and 
incorporate the inter-chain coupling afterwards. Each 
chain is then described by its own staggered magnetiza- 
tion field mi (i = 1,2) and its own nonlinear a model 
with Lagrangian density (||). After rescaling the fields as 
in Eq. ( pl| ) , we find that the inter-chain interaction takes 
the form C\ = hmi ■ m 2 , with h — 2Ks/a. In order to 
estimate the gap let us, as before, introduce Lagrange 
multiplier fields for the constraints m 2 = m\ = 1. The 
total Lagrangian density in imaginary time now reads 



I 2 

C E = - 22 [(dm,) 2 + (Tj(m^ - l/g)] + him • m 2 



(27) 



We have again set the characteristic velocity to 1. Wc 
should now integrate out mi and rri2, except that we are 
dealing here with a 'mass matrix' 



E = 



(Ti h 



(28) 



h a 2 

The eigenvalues of this mass matrix are 

A± = i {ai +a 2 ± y/(ai - o 2 ) 2 + 4/i 2 } (29) 

The effective potential V(a%, a 2 ) obtained by integrat- 
ing the fields mj^ is 

V(a 1 ,a 2 ) = ~(a 1 + a 2 )-~ / — — trln [lu 2 + k 2 + £ 
g 2 J 2ir 2n 

(30) 

The saddle-point equations dV/ do\ = and dV/ da 2 = 
may be expressed as 



- - — 1 A ' 



g 4tt v A7aT 

3 Q-i-o-2 A + 
°=16^A— A7 ln Al 



(31) 
(32) 



The solution to these equations is oj = a 2 = a and 
A± = a ± h, with 



a = a/ Co + h 2 



(33) 



where cto is the h = value. The square of the gap to 
the lowest excited state being A_, this gap, relative to its 
h = value, is 



= yVl + h 2 -h 



(34) 



where h = H/Aq. If wc substitute for A the saddle- 
point value and the known numerical value 0.41J, we 
find h « 23.8K/J. According to the relation (||), the 
gap decreases linearly as K/ J increases and then flattens 
towards zero. Because of the large numerical factor (~ 
23.8) between h and K/J, this drop is quite rapid. Note 
that Eq. (34) is obviously wrong when K/J is too large. 



A/A 




FIG. 2. Reduced gap A/Ao as a function of inter-chain 
coupling K/J, as obtained in the saddle-point approximation 
of the coupled non-linear sigma models. The upper curve 
illustrates Eq. (24), while the lower curve is a numerical solu- 
tion of Eq. (41), obtained in a different saddle-point approxi- 
mation. 

The saddle-point result (|34|) will receive corrections in 
1/N. We shall not get into these calculations. Instead, 
let us show how a slightly different result may be obtained 
by applying the saddle-point approximation in terms of 
the masses square A± instead of o\. 2 . Let us make a lin- 
ear transformation within the internal space spanned by 
the fields mi and m 2 such as to make the mass matrix 
E diagonal. We proceed to a functional change of inte- 
gration variables, replacing cr.; by A±, and m, by m±. 
The Jacobian associated with this change of variables is 
trivial in the latter case, but not in the former: 



det 



da. 



A+-A- 



A/CA+ - A_ + 2h)(X+ - A_ - 2h) 



(35) 



If we call this Jacobian J(x,t), the partition function for 
this systems reads 



Z = 



! J d\± dm± J(x, t) 



cxp —Se 



(36) 



In order to bring this expression into its usual form, we 
must bring the Jacobian within the exponential: 
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l[j(x,T)=vq>-W[\ ± } 



where 

W[\±] = 



dxdr^ln 



(A + -A^) 2 -4fe 2 
(A+-A_) 2 



(37) 



(38) 



Here we have introduced the lattice spacing a in the imag- 
inary time direction as well as in the spatial direction, in 
order to incorporate the functional Jacobian into a local 
actiono We may now proceed to the functional integra- 
tion of the fields m±, which is done exactly like in the 
ordinary nonlinear a model, except that a is replaced by 
A±. The effective potential in the large- N limit is then 



X+ + X- 



1 



2g 2a 2 
3 \ f dk dio 



In 



(A+-A_) 2 



Ah 2 



2 ^ I 2vr 2vr 



In 



(A+-A_) 2 
A,; 



k 2 + i 



(39) 



We now want to solve the saddle-point equation for this 
potential. To this end, let us write A± ~ a ± /3. The 
saddle-point equations are 



dV 
da 
dV 
dp 







= = - 



1 3 , 
- + — In 

g Air 

h 2 



A 2 



a 2 (3{(3 2 - h 2 ) 



3 lri a + P 

8ir a — [3 



(40a) 
(40b) 



Here again A is a momentum cutoff of the order of I /a. In 
the decoupled limit [h — 0) the solution to these equa- 
tions is (3 = and a = a Q = 2Aexp — 4ir/3g. If we 
now define reduced variables a = a/ao, f3 = /3/ao and 
h = h/ao, we can put the above saddle-point equations 
in the following form: 



P 2 



1 



f3 { p 2 -h 2 )ln^±i = r, 2 h 2 
a — p 



(41a) 
(41b) 



where rj 1 = 87r/3aoa 2 . The gap at finite inter-chain cou- 
pling is then A (AT) = A(0) v / a - 0. The equations (|l|) 
may be solved numerically. We need some input about 
the parameter rj: using the saddle-point value for ao, one 
finds rj ~ 14.1. A numerical solution of Eq. (|IITj) is illus- 
trated on Fig.(||), along with the relation (|34|). Both of 
these results are obtained in the saddle-point approxima- 
tion, but in terms of different variables. Of course, the 
change of variables defined in Eq. ( |29| ) should not make 
any difference if the exact value of the gap could be cal- 
culated; the 1/N expansions in terms of A± and CT1.2 are 
however different. 



IV. NUMERICAL RESULTS 

In this section we present some numerical results ob- 
tained for the spin-1 Heisenberg ladder (s = s = 1) and 



compare them with the analytical predictions worked out 
in the previous two sections of the paper. 

The system ([!]) may be studied numerically by essen- 
tially two techniques: exact diagonalization and quan- 
tum Monte-Carlo. Exact diagonalizations may further 
be refined by applying the density-matrix renormaliza- 
tion procedure. The disadvantage of exact diagonaliza- 
tions is of course that the amount of computer mem- 
ory and computer time needed grows exponentially with 
system size, whereas Monte-Carlo simulations are time- 
consuming but may easily be implemented on medium- 
size systems. 

Before presenting results, let us briefly describe the 
method. We used a projector Monte-Carlo-XPMC) tech- 
nique, adapted from the work of Takahashilij, with slight 
modifications. The method requires that the off-diagonal 
elements of the Hamiltonian be non positive. The Heisen- 
berg Hamiltonian ([!]) does not immediately fulfills this 
condition, but on a bipartite lattice such as the ladder, 
a unitary transformation H — > UHU^ 1 may be defined 
that brings the Hamiltonian into this form, with 



(42) 



Let us still denote the resulting Hamiltonian by H . 

The Hilbert space of the system may be described by 
the basis of eigenstates of {S*,S*}: let us write such a 
basis state as |{si})- The essence of the PMC method 
is the discrete (imaginary) time evolution of an ensemble 
of walkers. A walker is a state that belongs the the basis 
described above, changing from time to time according to 
stochastic transitions. Each walker \n) is assigned a non 
negative weight w n (r) which depends on time. The walk- 
ers are evolved by repeated application of the operator 
W = exp —eH, where e is some small time interval. This 
evolution operator is applied stochastically, meaning that 
upon application, a state |{si}) is changed to an other 
state |{sj}) with a probability proportional to the matrix 
element ({sJ}|W|{si}). At the same time, the weight w n 
is multiplied by the diagonal element ({sj}| W|{sj}). Of 
course, the exponentiation of H is impossible to achieve 
exactly without knowing the exact solution to the prob- 
lem. We therefore divide the Hamiltonian into different 
parts, each being the sum of mutually commuting two- 
spin Hamiltonians. For the ladder, one needs three such 
parts, as follows: 



Hi — J |^Sj • Sj + i + Sj • Sj_|_i 

i even 

H2 = J |^Sj • Si+i + Sj • Si+i 

i odd 



H 3 = Kj2Si-Si 



(43a) 
(43b) 
(43c) 



The discrete evolution operator is then replaced by W — 
W1W2W3, where Wi — exp— eHi. Since the different 
terms in Hi commute, we only need to exponentiate a 



G 



two-spin interaction, which can be done analytically. The 
smaller e is, the smaller the error made in this approxi- 
mation to the real exponential. 

The iteration starts with a random distribution of N 
walkers, all with equal weights (e.g. w n = 1). More pre- 
cisely, the walkers are chosen randomly within the sub- 
space of interest, i.e., for a fixed value of S^ ot . The num- 
ber N is typically 10 000. At each step the imaginary 
time is increased by e (typically e ~ 0.05 J) and the op- 
erator W is applied as described above. After a few iter- 
ations the distribution of weights becomes too sparse (its 
standard deviation becomes too large in comparison with 
the mean) and the weights have to be reconfigured by a 
procedure eliminating the ligth walkers (see Refs. [l6| and 
|l7| ) and normalizing the weights. This reconfiguration 
procedure is applied regularly, but should not affect the 
physical properties of the ensemble when done properly. 
After the system has evolved for some time r = n T e, the 
energy of the ensemble of walkers is evaluated by project- 
ing onto some arbitrary state For convenience, this 
state is taken to have an equal projection on each of the 
basis states |{sj}). The evolution then continues, an en- 
ergy measurement begin taken every n T iterations. The 
time interval n T should be short enough to allow as many 
measurements as possible, but long enough for successive 
energy measurements to be statistically uncorrelated. 



eliminated by projecting onto the uniform state 

The outcome of the simulation for the gap and the 
ground state energy per site as a function of K is pre- 
sented on Figs. 3 and 4. Notice that the gap exhibits 
a sharp drop until about K ~ 0.5, after which it rises 
slowly. We have not performed simulations on ladders 
longer than 24 rungs (48 sites), which is not enough to 
do a satisfactory finite-size analysis for most values of K. 
Nevertheless, we feel that the correct qualitative behavior 
as a function of K may be extracted. 
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FIG. 3. Results of a Quantum Monte Carlo simulation for 
the gap A divided by the K = value Ao. The data are 
shown for N — 4, 8, 16 and 24, where N is the length of the 
ladder. Ao was also obtained from exact diagonalizations, 
when available. 

In order to measure the spin excitation gap this way, 
one first finds the lowest energy level in the Sf ot = 
sector, and then in the S* ot = 1 sector. The difference 
should be equal to the gap A. Notice that the unitary 
transformation ( [42| ) has effectively shifted the momen- 
tum of excitations by tt, so that the first excited state 
of a Haldane system has now momentum zero and is not 



FIG. 4. Results of a Quantum Monte Carlo simulation for 
the ground state energy per site eo = Eo/N, with N = 4, 8, 
16 and 24. 



V. DISCUSSION 

In the first part of this work we mapped the Hamil- 
tonian of the Heisenberg ladder to the one-dimensional 
nonlinear sigma model, within the path integral formal- 
ism and the approximation of local AF order. From our 
knowledge of the nonlinear sigma model (with or without 
topological term), we concluded that the ladder has an 
excitation gap for all values of the inter-chain coupling K 
if the spins s and s of the two chains are both integers, or 
both half- integers (i.e. if s + s is an integer), whereas it 
has no gap if s+s is a half-integer. In the former case, the 
ladder is in a Haldane phase for all nonzero values of K, 
positive or negative. This agrees with the interpretation 
of Ref. 

However, the quantitative estimate (|24| ) of the gap A 
as a function of A, obtained via the saddle-point approx- 
imation of the nonlinear sigma model and illustrated on 
Fig. 1, is clearly wrong in both the K — > and K — > oo 
limits. For small K, it announces a uniform increase of 
the gap from a finite value at K = 0. In reality, the gap 
vanishes as K — > if s — s =ps, whereas it drops from 
its K — value if s = s = l.fc3 Thus, there is a quali- 
tative difference between half-integer and integer spin in 
the ladder when K is small, as we naturally expect. This 
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difference does not show up in the nonlinear sigma model 
representation of the ladder. 

For large values of K, Eq. ( pi] ) predicts a square-root 
behavior ~ \[K when K is large. In reality, the gap is 
approximately equal to K in that limit, i.e., the energy 
difference between a singlet and a triplet on a single rung, 
when J is neglected before K. In fact, we do not expect 
the parameters of the ladder non-linear sigma model of 
Sect. II to be correct when the ratio K/J is too far from 
unity. However, for moderate values of K (0.5 < K < 2), 
the gap of the spin-1 ladder illustrated in Fig. 3 shows 
a moderate increase whose slope is compatible with the 
estimate of Fig. 1. 

The behavior of the ladder for small inter-chain cou- 
pling (K <C J) is better represented by two interacting 
nonlinear sigma models. In the half-integer spin case, the 
presence of a topological term in each of the chain makes 
the analysis quite difficult; we hope to report on this in 
later work. However, the integer-spin case lends itself 
to a saddle-point evaluation of the gap which is much 
closer to reality. In the absence of variational principle 
or other guideline, it is difficult to say which of the two 
saddle-point methods illustrated on Fig. 2 offers a bet- 
ter estimate of the gap. A comparison with the Monte 
Carlo data of Fig. 3 tilts the balance in favor of the lower 
curve, which displays a sharper drop as K increases from 
zero. Recall that the saddle-point approximation leading 
to this curve neglects the fluctuations of the mass matrix 
E of Eq. (Eq), whereas the upper curve is obtained by ne- 
glecting the fluctuations of the Lagrange multipliers (7$ of 
Eq. (If). 
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* Choosing the same lattice spacing in the imaginary-time 
and space directions is somewhat arbitrary, but does the 
least violence to the Lorentz invariance of the theory. 
** If the numerical results of Ref. ^ for finite- length ladders 
were plotted in terms of A/Ao, as in Fig. 3, we would 
observe a rise of the curve as the length L increases, in 
contrast to Fig. 3. 
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